home *** CD-ROM | disk | FTP | other *** search
Text File | 1994-06-05 | 4.7 KB | 198 lines | [MATS/MATL] |
- echo off;
- % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
- % To accompany the text:
- % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
- % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
- % This free software is complements of the author.
-
- % Algorithm 4.2 (Polynomial Calculus).
- % Section 4.2, Introduction to Interpolation, Page 212
- echo on; clc; format short; hold off; clear
-
- % Proceed with an investigation of synthetic division for
-
- % evaluating Pn(x), P'n(x) and ∫Pn(x)dx where
-
- % Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)
-
- % It is assumed that the constant K of integration is K = 0.
-
- % Pn(x), P'n(x) and ∫Pn(x)dx are evaluated by invoking the Matlab
-
- % subroutine polyval(C,X) polyval(Cd,X) polyval(Ci,X), respectively.
-
- pause % Press any key to continue.
-
- clc;
- % Several Taylor polynomial coefficient lists are
-
- % stored in M-files named; zcos zsin ztan zexp
-
- % zacos zasin zatan zcosh zsinh zsqrt zlog
-
- % zsqrt4 zinv zemx2d2 zcosde zsinde zlogq
-
- % The user must determine the appropriate formula
-
- % for f'(x) in order to graph f'(x) and P'n(x).
-
- % The user must also determine the formula for ∫f(x)dx
-
- % in order to graph ∫f(x)dx and ∫Pn(x)dx.
-
- pause % Press any key continue.
-
- clc;
- % The algorithms are illustrated with the coefficients of the
-
- % Taylor polynomial of degree n = 11 for f(x) = ln(1+x).
-
- % Issue the command zlog to load the coefficients
-
- % into the array C. The function name is loaded
-
- % into the variable fun, the maximum degree is loaded into N.
-
- % For graphing over [a,b] the endpoints are in a and b.
-
- % Load the Taylor coefficients.
-
- [fun,dfun,ifun,x0,m,C,Ax] = zlog;
-
-
- pause % Press any key continue.
-
- clc;
- a = Ax(1,1); % You can change the endpoint a.
- b = Ax(1,2); % You can change the endpoint b.
- c = Ax(1,3);
- d = Ax(1,4);
- n = 11; % You can change the degree n.
- C = flipud(C);
- C = C(1:n+1);
- C = flipud(C);
-
- pause % Press any key to graph f(x) and P(x).
-
- clc; clg;
- h = (b-a)/200;
- X = a:h:b;
- x = X;
- Y = eval(fun);
- Z = polyval(C,X);
- axis([a b c d]);...
- plot(X,Y,'-g',X,Z,'-r');...
- hold on;...
- plot([a b],[0 0],'b',[0 0],[c d],'b');...
- xlabel('x');...
- ylabel('y');...
- Mx1 = [fun,' and P'];...
- Mx2 = [Mx1,num2str(n),'(x)'];...
- title(Mx2);...
- grid;...
- axis;...
- hold off;...
- shg; pause % Press any key to continue.
-
- Mx1 = 'The function is f(x) = ';
- Mx2 = 'The interval is ';
- Mx3 = 'Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
- Mx4 = 'The degree is n = ';
- Mx5 = ', and the coefficients list c is:';
- clc,echo off,diary output,...
- disp(''),disp([Mx1,fun]),...
- disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
- disp(Mx3),disp([Mx4,num2str(n),Mx5]),...
- for i=1:5:n+1, disp(C([i:min(i+4,n+1)])'); end,...
- diary off, echo on
-
-
-
-
-
- pause % Press any key to graph f`(x) and P`n(x).
-
- clc; clg;
- a = Ax(2,1);
- b = Ax(2,2);
- c = Ax(2,3);
- d = Ax(2,4);
- x = X;
- Y = eval(dfun);
- Dv = zeros(n,1);
- for j=1:n,
- Dv(j) = n-j+1; % Form vector [n,n-1,...,1]
- end
- Cd = C(1:n).*Dv; % Coefficients for Pn`(x).
- Z = polyval(Cd,X);
- axis([a b c d]);...
- plot(X,Y,'-g',X,Z,'-r');...
- hold on;...
- plot([a b],[0 0],'b',[0 0],[c d],'b');...
- xlabel('x');...
- ylabel('y');...
- Mx1 = [dfun,' and P`'];...
- Mx2 = [Mx1,num2str(n),'(x)'];...
- title(Mx2);...
- grid;...
- axis;...
- hold off;...
- shg; pause % Press any key to continue.
-
- Mx1 = 'The function is f`(x) = ';
- Mx2 = 'The interval is ';
- Mx3 = 'P`n(x) = d(1)x^n + d(2)x^(n-1) + ... + d(n)x + d(n+1)';
- Mx4 = 'The degree is n = ';
- Mx5 = ', and the coefficients list Cd is:';
- clc,echo off,diary output,...
- disp(''),disp([Mx1,dfun]),...
- disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
- disp(Mx3),disp([Mx4,num2str(n-1),Mx5]),...
- for i=1:5:n, disp(Cd([i:min(i+4,n)])'); end,
- diary off, echo on
-
-
-
-
- pause % Press any key to graph ∫f(x)dx and ∫Pn(x)dx.
-
- clc; clg;
- a = Ax(3,1);
- b = Ax(3,2);
- c = Ax(3,3);
- d = Ax(3,4);
- x = X;
- Y = eval(ifun);
- Iv = zeros(n+1,1);
- for j=1:n+1,
- Iv(j) = 1/(n-j+2); % Form vector [1/(n+1),...,1/2,1]
- end
- Ci = C(1:n+1).*Iv; % Coefficients for ∫Pn(x)dx.
- Ci = [Ci;0];
- Z = polyval(Ci,X);
- axis(Ax(3,:));...
- plot(X,Y,'-g',X,Z,'-r');...
- hold on;...
- plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
- xlabel('x');...
- ylabel('y');...
- Mx1 = [ifun,' and ∫P'];...
- Mx2 = [Mx1,num2str(n),'(x)dx'];...
- title(Mx2);...
- grid;...
- axis;...
- hold off;...
- shg; pause % Press any key to continue.
-
- Mx1 = 'The function is ∫f(x)dx = ';
- Mx2 = 'The interval is ';
- Mx3 = '∫Pn(x)dx = a(1)x^n + a(2)x^(n-1) + ... + a(n)x + 0';
- Mx4 = 'The degree is n = ';
- Mx5 = ', and the coefficients list Ci is:';
- clc,echo off,diary output,...
- disp(''),disp([Mx1,ifun]),...
- disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
- disp(Mx3),disp([Mx4,num2str(n+1),Mx5]),...
- for i=1:5:n+2, disp(Ci([i:min(i+4,n+2)])'); end,
- diary off, echo on
-